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I. Introduction 

There has been some recent work to develop two and three-dimensional alternating direction 
implicit (ADI) FDTD schemes [ 1 ]— [4] . These ADI schemes are based upon the original ADI concept 
developed by Peaceman and Rachford [5] and Douglas and Gunn [6], which is a popular solution 
method in Computational Fluid Dynamics (CFD). These ADI schemes work well and they require 
solution of a tridiagonal system of equations. A new approach proposed in this paper applies a 
LU/AF approximate factorization technique from CFD to Maxwell’s equations in flux conservative 
form for one space dimension. The result is a scheme that will retain its unconditional stability in 
three space dimensions, but does not require the solution of tridiagonal systems. The theory for this 
new algorithm is outlined in a one-dimensional context for clarity. An extension to two and three- 
dimensional cases is discussed in [7]. Results of Fourier analysis are discussed for both stability and 
dispersion/damping properties of the algorithm. Results are presented for a one-dimensional model 
problem, and the explicit FDTD algorithm is chosen as a convenient reference for comparison. 

II. One-Dimensional LU/AF Algorithm 


To develop the one-dimensional LU/AF algorithm. Maxwell’s equations for the one-dimensional 
TE mode are written in flux conservative form to include the electric and magnetic conduction 
currents as 


dq df_ 
dt dx 


-Pq 


( 1 ) 


where P is given by 


P = 


a /e 0 

0 <T*/p, 


( 2 ) 


and q = [E y , H Z ] T , f = [H z /e , E y //j] T and T denotes transpose. To develop the upwind LU/AF 
algorithm, the flux conservative form of Maxwell’s equations in (1) is recast in the following form 


dq -dq - 
dt + A d x ® 


where A is the Jacobian of / = Aq and is given by 


.4 = 


df 

dq 
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(3) 


(4) 


The matrix A has eigenvalues Ai, 2 = ±1 j s/Jdl corresponding to right and left propagating waves 
with speeds ±c. The eigenvalue matrix is given by 


A = 


\\ 0 

o a 2 


(5) 


The matrix A can be obtained from the eigenvalue matrix via a similarity transformation given by 
A = SAS 1 where the matrix S is composed of the eigenvectors of A The eigenvalue matrix 
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A can be split into two parts, one each for the right and left propagating waves and is given by 
A = A+ + A - . Using this relation, we can split A into two parts as _4 ± = SA^S -1 . The flux 
vector / is then split into two parts given by / = f + + f~ and the relation f + = A + q. For 
more details on this development, see reference [7]. This flux-vector-splitting method is similar to 
that developed by Steger and Warming [8] for the Euler equations governing inviscid fluid flow. To 
construct the LU/AF algorithm, time and space are discretized to t n = nAt, x , = iAx and the time 
derivative in (3) is approximated by a /i-weighted, O (At 2 ) difference equation given by 

dq „ (2/j + l)g" +1 - 4/jgf + (2/3 - l)#" 1 
dt ~ 2 At 

The spatial derivative is again replaced by a /(-weighted average between time level n -F 1 and n. 
This can be rewritten using the operator notation as 

g * P (A -f+ + A ijr) n+1 + (1 - P) (A ,Jt + A +fr) n (7) 

The A^ operators are now ( ) ( Ax 2 ) difference operators to be defined shortly. The parameter P can 
be used to construct a series of explicit and implicit schemes. For example, if p = 0, this results in a 
leapfrog scheme; p = 0.5 results in a Crank-Nicolson scheme and p = 1 results in an Euler implicit 
scheme. Using (6) and (7), the finite-difference equation for (3) is 

W + 1 -»**-' + „ (Aiff + At./D” + ' + 

(1 - ffl (Ai# + A +/f )” = -up ,-;■+• - (l - :j)pg (8) 

This can be rearranged as 

\l/(2At) + PaP + Pa (A^A+p) + A+ A ~ (■))] Aq" = -BP, (9) 

where A q? = g" +1 — g^ 1 - The residual, Rl^, is defined by 

mi = a jpg” + A. J; .. l g" + A+ A-q? - ‘^-A q?- 1 } (10) 

wither = 1/ ( 2/1 f- 1 ) . The difference operators in the residual are replaced by 0{Ax 2 ) backward and 
forward upwind difference operators on three-point one-sided stencils defined by 

AT ( (-) = (3(-)* - 4(-)i-i + (Ot- 2 ) /(2Ax) and A^(-) = ( 3(», + !(•), . , (•)/ • 2 ) /(2A,-;. 

Equation (9) is an 0(At 2 , Ax 2 ), unfactored, upwind scheme for electromagnetics. The LU/AF 
scheme is defined by factoring the left side of (9) into two operators, each designed for a forward 
and backward grid sweep as in the first order implementation. The LU/AF scheme is then given by 

[l/(2At)+paP + paAp i A+(-)]Aq* = -P" (11) 

\l/(2At) +paP + paA+A-{-)] Aq? = [//(2Ai) + paP] Aq* (12) 

III. Fourier Analysis 

A Fourier analysis shows that both the first and second-order upwind LU/AF algorithms are un- 
conditionally stable for p > 1/2, and that they contain both numerical dissipation for damping) 
and dispersion. The dissipation is present due to even order spatial derivatives in the truncation 
error which are a result of the upwind approximation. The complete Fourier analysis for dispersion 
and damping errors showed that the LU/AF algorithm has larger dispersion errors than the explicit 
FDTD method for v < 1. This is not particularly troublesome, because for i> < 1, explicit schemes 
are generally more efficient and are preferred. However, for v > 1, explicit schemes cannot be used 
and the LU/AF algorithm has its lowest dispersion error around v = 2. The results also showed that 
the lowest numerical dispersion is obtained when P = 1/2. The dispersive errors also increased with 
P and //. Since the LU/AF method uses windward differencing, numerical dissipation (or damping) 
is present in the solution. The dissipation is very small using a second-order LU/AF algorithm, and 
it decreases as the grid resolution is increased or as v is increased. 



IV. Boundary Conditions 


A. Outer Radiation Boundary Condition 

The characteristic based LU/AF method requires no extraneous boundary condition such as the 
Liao absorbing boundary condition [9] or the PML [10]. The only additional information required 
is information about waves that are entering the domain. Waves exiting the domain are naturally 
handled by the interior point algorithm. Therefore, the characteristic boundary conditions are im- 
plemented as follows: at grid point i = 0, equation (12) is used along with a specification of the 
incoming, right-going flux, Jq. At grid point i = imax, equation (11) is used along with a speci- 
fication of the incoming, left-going flux, ff max - The only additional information introduced at the 
boundary is nothing more than what is required by the physical system. In multidimensional prob- 
lems, the local coordinates at the outer boundaries are rotated to align with the direction of wave 
propagation defined by E x H. The characteristic equations are developed along this direction and 
are appropriately applied at the boundaries. This procedure was discussed and outlined by Shang 
[ 11 ]. 

B. Dielectric Surface Boundary Condition 

Since the LU/AF scheme follows the direction of information propagation (i.e. the characteristic ), 
at a material interface, the slope of the characteristic curve (i.e. the speed of light in the material) 
changes. Therefore, for the LU/AF method to be widely applicable, a careful treatment of material 
interfaces is required. With a material boundary in place, the right-going and left-going character- 
istics see a change in characteristic speeds, and therefore, a material interface condition needs to be 
implemented to correctly model the physics. This feature was recently developed and is outlined in 
[ 12 ]. 


V. Results 

The second-order LU/AF algorithm was tested by implementing equations (11) and (12) for inte- 
rior grid points away from absorbing boundaries and a first-order scheme at grid points next to the 
absorbing boundaries. Characteristic-based boundary conditions were used to terminate the com- 
putational domain and the incoming flux (f irnax ) at the right boundary was set to zero. The code 
was initialized by writing a time snapshot of a propagating pulse in the grid. A standard Gaussian 

_ t \ 2 

pulse of the form E y j = e 1 r ' was used as an excitation source. The main parameters of 
interest are the time resolution (i.e. number of time steps/period) and the grid resolution (i.e. number 
of cells/wavelength) of the highest frequency of interest in any given problem. These parameters 
will be designated as A 7 t and N x , respectively. The highest frequency of interest, f max , is calculated 
based upon N x and the largest cell size, and the time step for the implicit scheme is calculated based 
upon f max and the desired time resolution, A 7 t . For the FDTD method, the time step was calculated 
based upon the Courant stability condition using the smallest grid size. 


To illustrate the use of the LU/AF scheme for EM 
problems, propagation on a uniform mesh and reflec- 
tion from a lossy dielectric half-space were considered. 
The dielectric interface scheme discussed in Section IV- 
B was implemented and the results are compared with the 
explicit FDTD method. To solve the reflection problem, 
the problem space is a uniform grid with 2000 cells, a cell 
size of 1 cm, and it is filled with a lossy dielectric mate- 
rial from cell numbers 751-2000. The time step is 33.3 
ps and 0 = 1. The total electric field is recorded versus 
time at cell number 750. The incident field is obtained 
by running the code with free space only and recording 
the field versus time at the same location. The incident 
field is subtracted from the total field to give the scattered 
field. A point source located at the left boundary of the 
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Fig. 1. Reflection coefficient magnitude versus 
frequency for scattering from a lossy dielectric 
half-space on a uniform mesh with p = 1/2 
and v = 1. 




grid was used as the excitation source and the dielectric half-space material parameters are e = 4eo. 
ji = ji(, and <7 = 0.2 S/m. Figure 1 shows the reflection coefficient results and the agreement is 
excellent. 

For the propagation problem, the pulse was allowed to propagate for approximtely 50 meters and 
periodic boundary conditions were used along with v = 2. Figure 2 shows the error in electric field 
for LU/AF method. We see that the results are accurate to within about 1%. 

From additional tests, it was determined that the LU/AF 
method produces less accurate results for moderate con- 
ductivity values. This is due to the factorization error 
terms. As a is increased, the factorization errors become 
larger. In principle, this factorization error can be reduced 
through iterative error reduction [7], For perfect conduc- 
tors, we simply set q = 0. 

VI. Conclusions 

This paper has introduced an implicit LU/AF FDTD 
method for computational electromagnetics. A second- 
order accurate, LU/AF, characteristic based algorithm for 
electromagnetics has been implemented and tested on one- 
dimensional model problems for uniform and nonuniform 
grids. The one -dimensional model problem results show 
that accurate solutions for wave propagation can be ob- 
tained using Courant number greater than one and the 
method works well with lossy dielectric materials. The present results demonstrate potential ad- 
vantages in a one dimensional context, and this approach appears promising for development of 
stable, accurate and efficient implicit LU/AF schemes for complex two and three dimensional appli- 
cations. Extensions to two and three-dimensional applications have been outlined [14], but details 
of this work will the subject of future reports and articles. 
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Fig. 2. Error in electric field versus distance for 
free space propagation on a uniform mesh 
with ft ~ 1/2 andzy = 2. 




